#ifndef __PDF__
#define __PDF__

#include "tensor.h"
#include "mtrand.h"
#include <iostream>
#include <vector>

using namespace std;

/*! 
  \brief PDF functions


  this class contains a table of templates
  values are interpolated between bins.
*/
  
class PDF{
 public:
  PDF(){valid=false;}
  //!< PDF must be initialize externally, a trivial
  //!< constructor is presented but PDF will be invalid
  //!< until constructed by DataSet class or read from a file

  double operator()(const vector<double>& ntuple, double* err=NULL) const;
  //!< evaluate the PDF at a location
  //!< using linear interpolation
  //!< also obtain error if a pointer is specified

  double FullEvaluate(const vector<double>& ntuple, double* err, 
		      vector<double>& grad_output) const;
  //!< evaluate the PDF, error and grab the gradient
  //!< using linear interpolation

  bool Write(const char* output) const;
  //!< function to output a PDF

  bool Read(const char* input);
  //!< function to read a PDF from file

  bool IsValid() const
  {return valid;}
  //!< returns true if PDF is valid

  int Rank() const
  {return data.Rank();}

  friend ostream& operator<<(ostream&, const PDF&);
  //!< streaming operator for debugging

  friend class DataSet;
  //!< allow DataSet to read private members

  Tensor<double> data;
  //!< store PDF value at a multidimensional grid

  /*! enum for choosing which axes to project out
   */

  enum Projection{
    KEEP, INTEGRATE
  };

  PDF Project(const vector<Projection> input) const;
  //!< to get a new PDF by projecting a given dimension

  /*! enum for setting which is input and out in MC_generator
   */
  
  enum MC_Flag{
    INPUT, OUTPUT
  };

  /*! \brief Monte Carlo Event


    the MC_Generator generates MC_Events,
    MC_events store weights and errors for each event
    the erros ARE NOT INDEPENDENT, treat them as another 
    variable and compute them by averaging, not adding
    in quadrature
   */

  class MC_Event;

  /*! \brief Monte Carlo Integrator


    this is used to generate monte carlo events based
    on the given PDF. Slices can be taken by the 
    appropriate settings
   */

  class MC_Generator;

  MC_Generator Generator(const vector<MC_Flag>& flag_) const;
  //!< returns a MC_generator given some flags


  vector<double> one_over_width;
  //!< vector to store width info

  vector<double> width;
  //!< vector to store width info

  vector<double> min;
  //!< vector to store min values

  vector<double> max;
  //!< vector to store max values

  vector<int> nbins;
  //!< vector to store number of bins

 private:
  
  Tensor<double> data_bias;
  //!< PDF value bias

  Tensor<double> data_var;
  //!< PDF value variance

  Tensor<double> data_err;
  //!< PDF value error

  bool valid;
  //!< to keep track of whether not the PDF is valid

  double min_err;
  //!< the error given to zero PDF. This is done since Poisson "error"
  //!< come in at low N
  //!< or specially, an asymmetric poisson confidence interval should be used
  //!< we patch things up by demanding a minimum absolute error applied
  //!< whenever the pdf is computed with very small statistics

  void Initialize(int nvar, const double* min, const double* max,
		  const double* precision_=NULL, const double min_err_=0);
  //!< initialize the PDF
  //!< PDF could still invalid!

  //private inline functions
  inline double GetBinCenter(int index, int bin) const;
  //!< compute value centered at bin (a particular dimension)
  
  inline vector<double> GetBinCenter(const vector<int>& index) const;
  //!< compute value array centered at an vectorial index 
  //!< i.e. pt, m, t21...etc

  inline int GetBin_helper(double value, int bin) const;
  //!< return bin given value
};


//more detailed implementation of the MC_Generator class
class PDF::MC_Generator{
 public:
  
  MC_Generator(const vector<MC_Flag>& flag_, const PDF* pdf_);
  //!< default constructor, fix the flags and PDF
  
  vector<MC_Event> Generate(const vector<double>& input, int nevt=1000) const;
  //!< generate a number of Monte Carlo events

  vector<MC_Event> GenerateGrad(const vector<double>& input, int nevt=1000) const;
  //!< generate a number of Monte Carlo events

  vector<MC_Event> Generate(int nevt=1000) const
    {
      return Generate(vector<double>(), nevt);
    }
  //!< generate a number of Monte Carlo events (without input)

  vector<MC_Event> GenerateFast(const vector<double>& input, int nevt=1000) 
    const;
  //!< generate a number of Monte Carlo events, much faster but binned

  vector<MC_Event> GenerateFast(int nevt=1000) const
    {
      return Generate(vector<double>(), nevt);
    }
  //!< generate a number of Monte Carlo events (without input)

  vector<MC_Event> GenerateFull(const vector<double>& input) const;
  //!< this does not give MC events, but rather simply dump the whole 
  //!< PDF into MC_Events

			 
 private:
  vector<MC_Flag> flag;
  //!< settings for which variables to generate MC
  
  const PDF* pdf;
  //!< the corresponding PDF

  int input_size;
  //!< size of the input
  
  vector<int> index;
  //!< which index to generate Monte Carlo

  static MTRand_int32 irand;
  //!< random number generator
  
};


//more detailed implementation of the MC_Event class
class PDF::MC_Event{
 public:
  
  double err;
  //!< error for the event

  double weight;
  //!< weight for the event

  vector<double> extra;
  //!< extra variable

  vector<double> value;
  //!< value for the monte carlo

  MC_Event(const vector<double>& value_, 
	   double err_, double weight_):
  value(value_), err(err_), weight(weight_) {}
  //!< default constructor

  void SetExtra(const vector<double>& extra_)
  {
    extra=extra_;
  }
  
  vector<double>& GetExtra()
  {
    return extra;
  }

  double operator[](int i) const
  {
    assert(i>=0 && i<value.size());
    return value[i];
  }
  //!< access operator (checks for bounds)

  

};


inline double PDF::GetBinCenter(int index, int bin) const
{
  //DEBUG
  /*
  cout<<"GetBin small input : "<<index<<","<<bin<<endl;
  cout<<"GetBin small output : "<<(index+0.5)*(max[bin]-min[bin])<<","<<min[bin]<<endl;
  */

  return (index+0.5)*width[bin]+min[bin];
}

inline vector<double> PDF::GetBinCenter(const vector<int>& index) const
{
  assert(index.size()==data.rank);
  vector<double> result(data.rank);
  
  //grab central value of bin in each dimension
  for(unsigned int i=0; i<result.size(); i++)
    {
      result[i]=GetBinCenter(index[i], i);
    }

  return result;
}


inline int PDF::GetBin_helper(double value, int bin) const
{
  double temp_bin= (value-min[bin])*one_over_width[bin];
  //if value is in underflowed, go to the zeroth bin
  if(temp_bin<0) return 0;
  
  //if value is overflowed, go to the last bin
  if(temp_bin>=nbins[bin]) return nbins[bin]-1;
  return static_cast<int>(temp_bin);
}


#endif
